Predicting the current habitat refugia of Himalayan Musk deer (Moschus chrysogaster) across Nepal

Abstract Himalayan Musk deer, Moschus chrysogaster is widely distributed but one of the least studied species in Nepal. In this study, we compiled a total of 429 current presence points of direct observation of the species, pellets droppings, and hoofmarks based on field‐based surveys during 2018–2021 and periodic data held by the Department of National Park and Wildlife Conservation. We developed the species distribution model using an ensemble modeling approach. We used a combination of bioclimatic, anthropogenic, topographic, and vegetation‐related variables to predict the current suitable habitat for Himalayan Musk deer in Nepal. A total of 16 predictor variables were used for habitat suitability modeling after the multicollinearity test. The study shows that the 6973.76 km2 (5%) area of Nepal is highly suitable and 8387.11 km2 (6%) is moderately suitable for HMD. The distribution of HMD shows mainly by precipitation seasonality, precipitation of the warmest quarter, temperature ranges, distance to water bodies, anthropogenic variables, and land use and land cover change (LULC). The probability of occurrence is less in habitats with low forest cover. The response curves indicate that the probability of occurrence of HMD decreases with an increase in precipitation seasonality and remains constant with an increase in precipitation of the warmest quarter. Thus, the fortune of the species distribution will be limited by anthropogenic factors like poaching, hunting, habitat fragmentation and habitat degradation, and long‐term forces of climate change.

The distribution of musk deer is influenced by several factors including habitat, climate, and anthropogenic aspects (Singh, Gautam, et al., 2020;Singh, Mainali, et al., 2020).Studies have reported that climatic variables have greatly contributed to the distribution of musk deer in the Nepalese Himalayas (Lamsal et al., 2018;Singh, Gautam, et al., 2020;Singh, Mainali, et al., 2020).Among the climatic factors, precipitation was recognized as the most important factor for predicting the habitat suitability of the species, on the suitable habitat was found in the higher precipitation areas (Khadka et al., 2017).Similarly, the temperature has also a prominent role in the habitat suitability of musk deer (Lamsal et al., 2018).The other important habitat variable is vegetation which is responsible for determining the habitat suitability of this species (Nandy et al., 2020).However, this vegetation growth or availability is positively associated with precipitation (Tiwari et al., 2017).Regarding the habitat variable distance to water sources, the probability of occurrence of musk deer decreased with the increase in distance to water sources (Thapamagar et al., 2021).Similarly, musk deer was found in the gentle slopes up to 20 degrees as mentioned by earlier studies (Aryal et al., 2010;Neupane et al., 2021).Regarding the anthropogenic factors, the musk deer usually avoids human activities and livestock grazing sites in human-dominated landscapes (Thapamagar et al., 2021).
However, there is a seasonal and temporary nature of settlements in the high Himalayas of Nepal, so the species might overlap with the settlements, particularly during the winter season when people move to the lowlands to avoid extreme cold weather (Nandy et al., 2020).
In species distribution modeling (SDM) or ecological niche modeling, the potential distribution of a species is explained by the relationship between the species and the surrounding ecological and environmental factors (Beery et al., 2021;Peterson & Soberón, 2012;Thuiller et al., 2009).SDM can be used as a conservation planning approach for threatened species by determining the species distribution range and ecological niche (Adhikari et al., 2019).Due to the presence of large data and multifaceted associations between species and ecological variables, the scope of computer algorithms for ecological niche modeling, habitat modeling, predictive habitat distribution modeling, and range mapping such as SDM has increased to solve the problem of ecologists and statisticians (Beery et al., 2021).Besides, SDM helps to envisage the effects of climate change on species, which is very essential to achieve the conservation goals of being aware of the species distribution (Forester et al., 2013;Raymond et al., 2020).Discrepancies among different SDMs create challenges in determining the optimal model choice (Elith et al., 2011;Elith & Leathwick, 2009;Renner & Warton, 2013).This becomes especially evident when models are employed to forecast species distribution in distinct scenarios, such as projecting it into varied geographic regions (Thuiller, 2004;Thuiller et al., 2009).The ensemble modeling approach offers a viable solution to navigate through these complexities.Through the ensemble method in SDM, several modeling techniques are assembled to improve the projecting performance (Hao et al., 2020).The temperature has been found as a significant variable in shaping the distribution of several Himalayan species (Elsen et al., 2017;Koju, Bashyal, & Shah, 2021;Koju, Chalise, & Kyes, 2021).So, the wildlife of higher elevations or mountainous regions is more vulnerable to the impacts of climate change (Aryal et al., 2016;Elsen et al., 2020).Change in the vegetation composition and shift in the vegetation range have been documented from different regions of the Himalayas: i.e. west (Lamsal et al., 2018), East (Manish et al., 2016), and Central (Chhetri & Cairns, 2015) which are major consequences of climate change.This anticipated climate change will alter the climatic niche and shift the geographical ranges of several faunal species in the future.For example, about 30% of snow leopard (Panthera uncia) living space is anticipated to be lost in the entire Himalayan area by 2050, of which 40% could vanish from Nepal (Forrest et al., 2012).Likewise, Aryal et al. (2016) anticipated diminished habitat for snow leopards and blue sheep (Pseudois nayaur) in Nepal in the future climate.All of these confirmations recommended that climatic change drives the species to modify their geographic distribution in each locale, including the Himalayas.With concern to HMD, despite the anthropogenic activities such as habitat loss, habitat degradation, and poaching being the major factors leading to the population decline of HMD (Harris, 2016;Jnawali et al., 2011;Neupane et al., 2021), the species are additionally affected by the effect of climate change with a fluctuating level of results over space and time (Lamsal et al., 2018;Van Gils et al., 2016), Thus, the conservation of such threatened and crucial species of Himalaya region is highly essential in the scenario of the projected increase in climateinduced warming on those regions.This study aims to evaluate the appropriateness of the habitat for the endangered HMD species by examining various ecological and anthropogenic factors in their existing distribution areas.It is crucial to identify suitable habitats to enhance connectivity and ensure long-term conservation efforts.We hypothesized that HMD might exhibit a preference for habitats near water sources while avoiding areas that have significant anthropogenic influence.Additionally, we anticipate that a substantial portion of potentially suitable habitat exists outside of the currently protected areas.By investigating these aspects, we aim to gain a comprehensive understanding of HMD's habitat requirements and contribute to effective conservation strategies.

| Study area
Nepal is a mountainous country that extends over 147,516 km 2 in South Asia between the latitudes of 26°22′ to 30°27′ north and longitudes of 80°04′ to 88°12′ east.Because of the variability of climate and topography along a strong altitudinal gradient spanning from 60 to 8848 m above mean sea level, the country is endowed with extensive biodiversity (Bhattacharjee et al., 2017;Paudel et al., 2012).
There are three major physiographic regions in Nepal: (1) lowland (Terai and Siwalik), (2) mid-hills, and (3) high mountains (Shrestha & Aryal, 2011).The prevailing climate in the country is characterized by dry winters and hot summers (Karki et al., 2016).The average annual precipitation is 1768 mm and the annual mean temperature is 18°C (Shrestha et al., 2000).The high mountains, which are the preferred distribution range of HMD, cover 24% of the country's total geographical area and comprise two-thirds of the country's PAs (Shrestha et al., 2010).

| Presence data of HMD
We obtained the occurrence locations of HMD mainly from the field-based surveys during 2018-2021 and Periodic data held by the Department of National Park and Wildlife Conservation of Nepal between 2018 and 2021.For field observations, we relied on direct sightings of the species, camera trapping as well as indirect sign sightings-pellets droppings and hoofmarks for determining the presence of HMD (Figure 1).Its pellets were differentiated from other sympatric deer species by their shape and size (Neupane et al., 2021).Within our pre-defined study period, a total of 429 existing presence points of HMD were compiled.
We used the spatial resolution of the environmental variables employed in this modeling was 1 km.Similarly, the SpThin package in R spatially attenuates the occurrence dataset ensuring that no two locations were inside a grid of 1 × 1 km (Aiello-Lammens et al., 2015).Thus, only one presence point in each grid cell was used to minimize spatial autocorrelation and avert inflated measures of accuracy (Veloz, 2009).Spatial filtering also helps to improve model prediction performance by reducing the effects of sample bias (Boria et al., 2014).Following filtering procedures, a set of 346 spatially independent HMD presence locations were retained and used for modeling.

| Environmental variables
We used a combination of bioclimatic, anthropogenic, topographic, and vegetation-related variables to predict the current suitable habitat for HMD in Nepal.Given that variable selection is regarded as a critical phase in SDM, an effort to incorporate important predictor variables (Araújo & Guisan, 2006) was used.
Initially, we identified a set of 33 variables (Table 1) based on literature, which suggested those variables were important for the habitat suitability of HMD.We performed the multicollinearity test for the selected environmental variables and avoided the environmental variables with correlation coefficients >0.8 and variance inflation factor (VIF) >5, which helped to prevent our model overfitting.Finally, 16 variables were retained as predictor variables in habitat suitability modeling for HMD as suggested by Zeng et al. (2016).

| Bioclimatic variables
For spatial modeling, bioclimatic variables are widely used given that these variables are ecologically important and characterize annual trends, seasonality, and temperature and precipitation extremes (Hijmans, 2012;Hijmans et al., 2005).WorldClim-Global Climate Data (www.world clim.org/ bioclim) was used to retrieve 19 bioclimatic variables (Fick & Hijmans, 2017).These data were retrieved in a grid format with a 1 km spatial resolution.

| Topographic variables
Topographic variables such as elevation, slope, aspect, and distance to water sources have governed the habitat suitability of megaherbivores (Ghimire et al., 2019;Sharma et al., 2020).In our study, elevation, aspect, and slope data were generated using ArcMap 10.8.1 (ESRI, 2020) with a 1 km spatial resolution Digital Elevation Model (DEM) acquired from the United States Geological Survey database (USGS, 2022).Shapefiles including water source information were downloaded from the Geobabrik website (GEOFABRIK, 2022) and transformed into a distance raster file using ArcMap10.8.1 (ESRI, 2020).

| Vegetation related variables
One of the most important elements determining the distribution of herbivores like HMD is vegetation-related variables (Gandiwa, 2014;Perea et al., 2015).The study collected four variables namely forest cover, minimum EVI, mean EVI, and maximum EVI.Forest cover was downloaded from Earth engine partner Appspot (Hansen et al., 2013).EVI time-series data was obtained from Moderate Resolution Imaging Spectroradiometer (MODIS) (USGS, 2022).In the TIMESAT algorithm (Jönsson & Eklundh, 2004), the Savitzky-Golay filter was employed to smooth the data.

| Predicting the distribution of HMD
In the first step, a multicollinearity test was conducted among the 33 environmental variables.Variables displaying a correlation coefficient >.7 and a variance inflation factor >5 were excluded to mitigate the multicollinearity effect (Dormann et al., 2013).As a result, 16 predictor variables remained and were utilized for habitat suitability modeling (Table 1).We followed the overview, data, model, assessment, and prediction (ODMAP) method suggested by Zurell et al. (2020) to create habitat suitability models for HMD in Nepal.
The utilization of ensemble maps in recent species distribution modeling (SDM) exercises has garnered considerable attention due to their demonstrated higher predicted accuracy (Hao et al., 2020).These maps are formed by merging multiple models constructed through diverse modeling approaches (Hao et al., 2019).
As a result, the habitat suitability model for HMD in Nepal was constructed using an ensemble modeling approach.The ensemble model was created in R (R Development Core Team, 2020) using the BIOMOD2 package (Thuiller et al., 2020) et al., 2008) and true skill statistics (TSS) (Allouche et al., 2006), are widely used measures to evaluate predictive performance (Thuiller et al., 2009).Despite being widely used as a measure for model evaluation, AUC is criticized for its shortcomings (Lobo et al., 2008).
Hence, the predictive performance of our model was evaluated using TSS criteria ranging between −1 and +1.To construct an ensemble model through a weighted mean strategy, all models having a TSS value >0.6 (Marmion et al., 2009) were selected.Three models (GBM, MaxEnt, and RF) have TSS value than 0.6 and hence we selected them to develop a weighted mean ensemble approach.
Secondly, we also analyzed the data with Maxent only model, to observe whether the parameter tuned maxent only model performs better compared to the ensemble model approach.The ENMeval package (Kass et al., 2021) in the R programming language was employed to optimize the MaxEnt model.In this study, a comprehensive evaluation was conducted on a set of 48 models.et al., 2022;Steen et al., 2019).Following the process of model optimization, the model that was selected as the best fit model had the characteristic class (Feature Class) LH, a value of RM equal to 1, and a delta AIC of 0. Following the adjustment of these parameter configurations, the maximum number of iterations was established at 1000 with 10,000 background points.Additionally, 70% of the presence points were allocated for training the model, while the remaining data was reserved for testing purposes (Barbet-Massin et al., 2012).
However, the accuracy of this Maxent only model (AUC-0.90,TSS-0.82) was less than the model accuracy obtained from Ensemble model (AUC-0.98,TSS-0.966).Therefore, we decided to keep the Ensemble model prediction to generate our results and discussion.

| Predicting the current suitable habitat of HMD
The current habitat suitability map generated through an ensemble modeling approach based on bioclimatic, topographic, vegetation, and anthropogenic variables indicated that 6973.76 km 2 (5%) area of Nepal is highly suitable and 8387.11km 2 (6%) is moderately suitable for HMD (Figure 2).
Within the overall suitable habitat of HMD, approximately 51.4% (7895.88 km 2 ) is located within protected areas of Nepal, while the F I G U R E 2 Present suitable habitat area for HMD (represented by green colors) across Nepal.
remaining 48.6% exists outside the protected area networks.Among the protected areas of Nepal, Khaptad National Park has the most suitable habitat (54.37%) followed by Dhorpatan Hunting Reserve (51.12%), while the Sagarmatha NP and Shey-Phoksundo NP has more suitable habitat in the buffer zones than inside the core area (Table 2).According to the administrative divisions, Karnali Province (3832.88km 2 ) has the most suitable area outside protected area networks followed by Sudurpaschim province (1184.6 km 2 ) (Figure 2; Table 3).

| Contributions of variables to build the model
Among the 16 predictive environmental variables used to predict current suitable habitat for HMD, climatic variables had greater impacts that described more than 90% of the model performance.
This model depicts Precipitation Seasonality (Bio15; 67%), followed by Precipitation of the Warmest Quarter (Bio18; 14%) as the major influencer for HMD distribution.Besides, the distribution of suitable habitats in the model was also influenced by the temperature ranges (both annual and diurnal) and distance to water bodies.And, among the anthropogenic variables, land use land cover change (LULC) had a major persuading factor in the model performance (Table 4).
The response curve of the models indicates that the probability of occurrence of HMD decreases with an increase in Precipitation Seasonality (Bio15) and remains constant with an increase in Precipitation of the Warmest Quarter (Bio18) (Figure 3).Similarly, the probability of occurrence of HMD in habitats far from water sources decreases continuously under Maxent while remaining constant after a certain distance under GBM and RF model.The probability of occurrence of HMD is less in habitats with low forest cover but the probability of occurrence remains constant with further increase in forest cover.

| Potentially suitable habitat for HMD
Our results indicate that suitable habitat for HMD is distributed throughout the mountainous regions of Nepal but is not thoroughly continuous, which aligns with previous studies conducted in the Nepalese Himalayas (Green, 1986).HMD has been recorded from all the Himalayan protected areas of Nepal (Aryal & Subedi, 2011).
Our model showed that 51.78% (8295 km 2 ) of the suitable habitat of musk deer is protected by the existing network of protected areas, which is in line with Lamsal et al. (2018), but contradicts with the prediction of Aryal and Subedi (2011), who suggested that only 19.26% (5815.08 km 2 ) of potential habitat lies within protected areas.However, the total potential habitat of HMD predicted by Aryal and Subedi (2011) is much larger than the 16,020 km 2 suitable habitat estimated in our study.Our model also revealed that 48.22% of suitable habitat lies outside protected areas, and the suitable habitat is not in a range throughout the mountainous regions which sustain the previous studies (Khadka et al., 2017;Lamsal et al., 2018;Singh, Gautam, et al., 2020;Singh, Mainali, et al., 2020).
Thus, we evaluated that identifying new areas outside PA networks is important for the conservation of HMD.Habitat distribution models provide bases to choose these areas of probable distribution of the species.Although HMD is one of the charismatic species of the Himalayas and is protected species under DNPWC Act, its studies and conservation actions are limited within the protected area.
We found that almost 50% of the probable habitat of HMD lies outside the protected area networks, thus it is of immense need to expand the focus of conservation actions beyond the Pas too.Similarly, the suitable habitat outside PAs is larger in the western landscape compared to the central-eastern complex.HMD is one of the most poached species (Ilyas, 2015;Subedi et al., 2012).Without prompt conservation actions and efforts beyond PAs networks, their survival will be increasingly confined to protected areas alone in the foreseeable future.Additionally, future climatically suitable habitat is predicted to be more in the western landscape of Nepal (Khadka et al., 2017).Although a good share of habitat lies within DHR and ANCA, research works are mere in these PAs compared to others.
Therefore, the areas in the western landscape of Nepal require serious attention for the survival and conservation of HMD.

| Distribution of HMD
Various factors such as climate, habitat, and anthropogenic variables influence the distribution of musk deer (Jiang et al., 2019;Singh, Gautam, et al., 2020;Singh, Mainali, et al., 2020;Wangdi et al., 2019).Our study's model predicted that climatic variables have a greater impact on the distribution of HMD in Nepalese Himalaya compared to habitat variables, which aligns with previous studies by Singh, Gautam, et al. (2020), Singh, Mainali, et al. (2020), and Lamsal et al. (2018).Precipitation is the most significant factor in determining the habitat suitability of HMD in Nepal Himalayas, which supports the findings of Singh, Gautam, et al. (2020) and Singh, Mainali, et al. (2020) but contradicts with findings of Lamsal et al. (2018) and Khadka et al. (2017).Specifically, our model predicted that HMD suitable habitat increases with an increase in precipitation during the warmest quarter, which is in line with research on other musk deer species (Singh, Gautam, et al., 2020;Singh, Mainali, et al., 2020).The model also suggested that the habitat suitability is influenced by precipitation in the dry winter and warm summer seasons, with a higher probability of occurrence of HMD in areas with higher seasonal precipitation, respectively (Figure 3).
In Nepal, the precipitation in summer is governed by monsoon, which brings heavy rainfall in the eastern part of the country and less rainfall in western regions, and the winter is governed by western disturbances which bring precipitation in the form of snow from the western region (Kansakar et al., 2004;Talchabhadel et al., 2018).The growth of vegetation in high mountains is affected by soil moisture, which in turn is influenced by temperature and precipitation patterns in the area (Paudel & Andersen, 2013;Regmi et al., 2020).The suitability of the HMD habitat is closely related to vegetation (Nandy et al., 2020).
The variability of precipitation seasonality across the region plays a significant role in determining the habitat suitability of HMD, as predicted by the model and reported in earlier studies (Khadka et al., 2017;Lamsal et al., 2018).As a result, the model projected a larger area of highly suitable habitats in the centraleastern landscape compared to the western landscape (as shown in Figure 2 and Table 3).Vegetation growth in the trans-Himalayan range was found to have a positive relationship with precipitation only (Ale et al., 2018;Tiwari et al., 2017).Studies have also reported that precipitation on the leeward sides is very low (Talchabhadel et al., 2018), which might be one of the reasons why the model predicted most of the suitable habitats in the gullies of river valleys.Overall, Nepal's precipitation seasonality is important in determining the distribution of HMD.Besides precipitation, the model diagnosed a strong response with temperature variables in determining the probability of occurrence in Nepalese Himalaya as demonstrated by mean diurnal range (BIO 2) and isothermality (BIO 3) being the third-most influential variables to the model.This relationship is consistent with previous studies (Khadka et al., 2017;Lamsal et al., 2018).
The likelihood of HMD occurrence was also influenced by topographic variables.Firstly, the proximity to water sources negatively impacts the probability of detecting HMD, with the likelihood of occurrence decreasing as the distance from water sources increases.
This finding is consistent with previous research that highlights the importance of water availability in determining musk deer distribution (Singh et al., 2018;Thapamagar et al., 2021).Secondly, our model also predicted that HMD prefers gentle slopes, which aligns with earlier studies that report an increase in the likelihood of musk deer occurrence with slopes up to 20°, followed by a slight decrease in probability beyond that point (Aryal et al., 2010;Neupane et al., 2021;Thapamagar et al., 2018).Figure 3

F I G U R E 3
Response curves indicating the effects of different environmental variables on habitat suitability of HMD.
be due to the seasonal and temporary nature of settlements in the high Himalayas in Nepal, where musk deer habitat might overlap with these settlements (Nandy et al., 2020).However, the model depicted that the probability of musk deer occurrence decreases far away from foot trails (Figure 3).

| Conservation implications
In

ACK N OWLED G M ENTS
We would like to thank the Department of National Parks and Wildlife Conservation (DNPWC) of Nepal for data support.We are grateful to the field assistants and local people who assisted us throughout the fieldwork.Finally, we would like to acknowledge the contributions of the officials of the protected areas (habitat ranges of HMD) for their cooperation and support required for this study.

CO N FLI C T O F I NTER E S T S TATEM ENT
None declared.

O PEN R E S E A RCH BA D G E S
This article has earned Open Data, Open Materials and Preregistered Research Design badges.Data, materials and the preregistered design and analysis plan are available at https:// doi.org/ 10. 5061/ dryad.sj3tx 969j.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data associated with this manuscript are available at https:// d at a d r y a d .o r g / s t as h/ s h a r e/ l F T Z j u S Kwc a L l a G KQ k m I Aot r Q 1caH62_ HZxgZ J1_ 60M.

F
I G U R E 1 Map of study area showing the protected area types within Nepal and occurrence distribution throughout the study area.
Current HMD suitable habitat within protected areas of Nepal.Overall suitable area by province, along with the total suitable area that falls outside the protected area network.